% SpeedyGA is a vectorized implementation of a Simple Genetic Algorithm in Matlab
% version 1.2.1
% Copyright (C) 2007,2008  Keki Burjorjee

%  Licensed under the Apache License, Version 2.0 (the "License"); you may
%  not use this file except in compliance with the License. You may obtain 
%  a copy of the License at  
%
%  http://www.apache.org/licenses/LICENSE-2.0 
%
%  Unless required by applicable law or agreed to in writing, software 
%  distributed under the License is distributed on an "AS IS" BASIS, 
%  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 
%  See the License for the specific language governing permissions and 
%  limitations under the License. 

% Created and tested under Matlab 7 (R14). 


% If a parameter is passed to the speedyGA function, it is used to set the
% seed of the random number generator. Otherwise, the seed is made to be
% dependent upon the system clock. 
% Sample usage: speedyGA(), speedyGA(27463)

function [eliteIndiv eliteFitness]=speedyGA(varargin)
len=640;                 % The length of bitstrings
popSize=500;             % The size of the population
maxGens=500;             % The maximum number of generations allowed for a run
probCrossover=1;         % The probability of crossing over. 
                         % (Note: 1-probCrossover is the probability of
                         % reproduction without crossover)
probMutation=1/popSize;  % The mutation probability (per bit)
sigmaScalingFlag=1;      % Sigma Scaling is described on pg 168 of M. Mitchell's
                         % GA book. It often improves GA performance.
sigmaScalingCoeff=1;     % Higher values => less fitness pressure 

crossoverType=2;         % 0 => no crossover
                         % 1 => 1pt crossover
                         % 2 => uniform crossover

visualizationFlag=1;     % 0 => don't visualize bit frequencies
                         % 1 => visualize bit frequencies

quietFlag=1;             % 0 => display details of each generation
                         % 1 => run quietly


%pregenerate crossover and mutation masks
maskReposFactor=5;
crossmasks=[];
if crossoverType==0
    crossmasks=false(popSize,len);
elseif crossoverType==2
    crossmasks=rand(popSize,(len+1)*maskReposFactor)<0.5;
end
mutmasks=rand(popSize,(len+1)*maskReposFactor)<probMutation;

% preallocate vectors for recording the average and maximum fitness in each
% generation
avgFitnessHist=zeros(1,maxGens+1);
maxFitnessHist=zeros(1,maxGens+1);

% If a parameter was passed to the vectorGA function, it is used to set the
% seed of the random number generator. Otherwise, the seed is made to be
% dependent upon the system clock
if nargin==1
    rand('state',varargin{:});
else
    rand('state',sum(100*clock));
end
    
eliteIndiv=[];
eliteFitness=-realmax;


% the population is a popSize by len matrix of randomly generated boolean
% values
pop=rand(popSize,len)<.5;

for gen=0:maxGens
        
     fitnessVals=oneMax(pop);
    
    [maxFitnessHist(gen+1), maxIndex]=max(fitnessVals);    
    avgFitnessHist(gen+1)=mean(fitnessVals);
    if eliteFitness<maxFitnessHist(gen+1)
        eliteFitness=maxFitnessHist(gen+1);
        eliteIndiv=pop(maxIndex,:);
    end
    
    % display the generation number, the average Fitness of the population,
    % the maximum fitness of an individual in the population, and an 
    % individual with the maximum fitness
    if ~quietFlag
        display(['gen=' num2str(gen,'%.3d') '   avgFitness=' ...
            num2str(avgFitnessHist(gen+1),'%3.3f') '   maxFitness=' ...
            num2str(maxFitnessHist(gen+1),'%.3d')]);
    end
    % Conditionally perform bit-frequency visualization
    if visualizationFlag
        figure(1)
        bitFreqs=sum(pop)/popSize;
        plot(1:len,bitFreqs, '.');
        axis([0 len 0 1]);
        title(['Generation = ' num2str(gen) ', Average Fitness = ' sprintf('%0.3f', avgFitnessHist(gen+1))]);
        ylabel('Frequency of the Bit 1');
        xlabel('Locus');
        drawnow;
    end    

    % if sigma scaling is on then perform it
    if sigmaScalingFlag
        sigma=std(fitnessVals);
        if sigma~=0;
            fitnessVals=1+(fitnessVals-mean(fitnessVals))/...
            (sigmaScalingCoeff*sigma);
            fitnessVals(fitnessVals<=0)=realmin;
        else
            fitnessVals=1;
        end
    end        
    
    
    % Normalize the fitness values and then create an array with the 
    % cumulative normalized fitness values (the last value in this array
    % will be 1)
    cumNormFitnessVals=cumsum(fitnessVals/sum(fitnessVals));
    
    % Preallocate an array for the post-selection population 
    newPop=zeros(popSize,len);
    % Perform fitness proportional selection using Stochastic Universal
    % Sampling
    SUSMarkers=rand(1,1)+[1:popSize]/popSize;
    gt1=SUSMarkers>1;
    SUSMarkers(gt1)=SUSMarkers(gt1)-1;
    [temp newPopIndices]=histc(SUSMarkers,[0 cumNormFitnessVals]);
    pop=pop(newPopIndices,:);

    %%%%%%%%%%%%%%% Perform uniform recombination and reproduction    
    % generate pairs of mating indices
    matingIndices=ceil(rand(popSize,2)*popSize);
    % create a new population using the first mate from each mating pair
    newPop=pop(matingIndices(:,1),:);
    % create a new population using the second mate from each mating pair
    secondMate(:,:)=pop(matingIndices(:,2),:);
    % create crossover masks for the entire population
    if crossoverType==0
        masks=crossmasks;
    elseif crossoverType==1
        masks=zeros(popSize, len);
        temp=ceil(rand(popSize,1)*(len-1));
        for i=1:popSize
            masks(i,1:temp(i))=1;
        end
    else    
        temp=floor(rand*len*(maskReposFactor-1));
        masks=crossmasks(:,temp+1:temp+len);
    end
    % generate an array of booleans indicating which first mates to leave 
    % untouched
    reprodIndices=rand(popSize,1)<1-probCrossover;
    masks(reprodIndices,:)=false;
    % implement uniform recombination with reproduction
    newPop(logical(masks))=secondMate(logical(masks));
    % implement mutation
    temp=floor(rand*len*(maskReposFactor-1));
    masks=mutmasks(:,temp+1:temp+len);
    pop=xor(newPop,masks);    
end
if ~quietFlag
    figure(2)
    %set(gcf,'Color','w');
    hold off
    plot([0:maxGens],avgFitnessHist,'k-');
    hold on
    plot([0:maxGens],maxFitnessHist,'c-');
    title('Maximum and Average Fitness Per Generation')
    xlabel('Generations')
    ylabel('Fitness')
end

% onemax
function fitness=oneMax(pop)
fitness=sum(pop,2)';    

% The royal roads function. The chromosome length (i.e. len) 
% should be a multiple of 8
function fitness=royalRoads(pop)
[popSize len]=size(pop);
fitness=zeros(popSize,1);
for i=1:8:len
    temp=sum(pop(:,i:i+7),2);
    temp=double(temp==8);
    fitness=fitness+temp;
end
fitness=fitness';
    





